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ABSTRACT 

Formation rates of compact-object binaries are often derived from population synthesis calculations. 
However, such calculations depend sensitively on a relatively large number of model input parameters. 
Given considerable uncertainty in those model parameters, the predicted inspiral rates for double 
compact objects relevant to gravitational-wave interferometric detectors have been shown to be are 
uncertain by several orders of magnitude. Typically, inspiral rates are estimated for only a small set of 
models with a remarkably poor coverage of the highly multi-dimensional parameter space (primarily 
because of limited computer resources) . Here, using as an example seven population-synthesis model 
parameters, we show that it is possible to derive fits of double-compact-object inspiral rates dependent 
simultaneously on all seven parameters. We find these fits to be accurate to 50% for binary black 
holes and to 40% for binary neutron stars. The availability of such fits implies that (i) depending on 
the problem of interest, it is not necessary to complete large numbers of computationally demanding 
population synthesis calculations; and (ii) for the first time, the sufficient exploration of the relevant 
phase space and the assessment of the uncertainties involved is not limited by computational resources 
and becomes feasible. 
Subject headings: Stars: Binaries: Close 



1. INTRODUCTION 

Over the last decade the question of the Galac- 
tic inspiral rate of binaries with two compact objects 
\ (neutron stars NS or black holes BH) has attracted 
attention primarily because of the development and 
planning of gravitational-wave interferometric detectors 
both on the Earth and in space (e.g., LIGO, VIRGO, 
GEO600 and LISA). Such rate estimates have been 
widely used in the assessment of gravitational inspi- 
ral detcctability, given assumed instrument sensitivities 
{Cutler & Thornc 2002). 

A number of different groups have calculated 
inspiral rates using population synthesis calcula- 
tions, most com m only with Monte Carlo methods 
jFrver et al. 1998'; 'Portegies-Zwart & Yungelson 199^ 
Bethc & Brown 1998; Frver Wooslcv & Hartmann 1999; 
Belczvnski Kaloger a fc Bulik 2002t 

Voss fc Tauris uch studies consider the 

complete formation history of double compact objects 
through long sequences of binary evolution phases. 
However, our current understanding of single and binary 
star evolution is incomplete. Therefore uncertainties 
are often parameterized and some of these parameter- 
izations are better constrained (usually empirically) 
than others. Depending on the binary type of interest, 
rate predictions are sensitive to different parameters 
and often degeneracies have been shown to exist. It is 
clear that a comprehensive a priori estimate of inspiral 
rates requires estimates of model uncertainties as well. 
However, current binary population synthesis codes 
are computationally demanding (depending on their 
level of sophistication) and covering a large part of the 
multi-dimensional parameter space has been proven 
very challenging. As a result, earlier studies have 
estimated model uncertainties by varying only one or 
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two parameters at a time, and considered at best a 
couple of dozen of models. 

In this Letter we discuss ways to reduce the compu- 
tational cost of population synthesis calculations (used 
here to focus on the formation rates of double com- 
pact objects). As a result, significantly more extended 
parameter-space searches become feasible. We make use 
of the basic concept of genetic algorithms, which are ideal 
for the exploration of highly multi-dimensional spaces. 
We validate the suggested ways and use them to exam- 
ine whether the derivation of fits for the inspiral rates as 
a function of a large number of free parameters is pos- 
sible. We find that, while strong correlations between 
model parameters are apparent, it is possible to fit inspi- 
ral rates to an acceptable accuracy (50% in this study, 
but possibly smaller if the original simulations used to de- 
rive the fits are more accurate). In what follows we first 
describe briefiy the population synthesis code we use in 
our analysis. We then discuss the basic concepts of ge- 
netic algorithms and describe how they can be used to 
increase the computational efficiency of synthesis codes. 
We conclude with the derivation of fits for NS-NS and 
BH-BH rates and discuss the implications of the avail- 
ability of such fits for deriving constraints on BH-BH 
rates in the future. 

2. STARTRACK POPULATION SYNTHESIS CODE 

To generate and evolve stellar populations until dou- 
ble compact-object formation occurs, we use the Star- 
Track code first developed by Belczynski, Kalogera, and 
Bulik (2002) [hereafter BKB] and recently significantly 
updated and tested as described in detail in Belczynski 
et al. 2004. 

Here we briefiy summarize the main parts of the com- 
putations relevant to this study. We generate a large 
number N of binaries specified without loss of generality 
by the mass mi of the primary (in solar masses Mq)] 
the mass ratio q = TO2/toi < 1 between the primary and 
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secondary objects; the seniiniajor axis A (in solar radii); 
and the orbital eccentricity e. We draw binaries from the 
distributions given in BKB Eqs. (2), (4), (5), and 

where Qc ~ 0.2 and r G [0, 3] is a model parameter; for a 
discussion see Kalogera & Webbink (1998).^ We adopt a 
uniform star formation rat e appropriate for our Galaxy 
l|Rana 199H iGilmore 200^ : each binary is assigned a 
birth time from a uniform distribution between the for- 
mation of our galaxy {t — 0) and the present (assumed 
at t — T — lOGyr). We evolve each binary until either 
it ceases to be a binary, or until the present, whichever 
comes first. 

During this evolution of N massive binaries, some num- 
ber n of a certain type of events occurs, corresponding to 
a rate of r = n/T in the simulation. To obtain the time- 
averaged event rate TZ = J {dn/dt)dt for the Galaxy, 
we scale this computational rate up by a scale factor s 

TZ = sx n/T (2) 

where s is a ratio of the number of stellar systems we 
have effectively simulated to the number of stellar sys- 
tems of the same type in the target system (here, the 
Milky Way). The appendix describes in detail how we 
determine this absolute normalization scale factor. 

Parameter Phase Space: For the study presented here 
we have chosen to vary seven (7) model parameters in the 
synthesis calculations. The choice is strongly guided by 
our past experience with double-compact-object popula- 
tion synthesis (BKB) and represent the model param- 
eters for which strong dependence has been confirmed. 
The list of these seven parameters is as follows (for more 
details see BKB): the fraction of transferred mass that 
is lost from the binary in phases of non-conservative 
mass transfer fa G [0, 1]; the common envelope efficiency 
coupled with the uncertainty in the donor star central 
concentration ax X G [0,1]; three parameters describ- 
ing the locations and relative weight of two Mawellians 
in a bimodal NS kick distribution {vi S [0, 200]km/s, 
V2 e [200, 1000]km/s, s G [0,1]); a parameter describ- 
ing the stellar wind strength w € [0, 1] relative to the 
reference-model assumptions; and the power-law expo- 
nent r g [0, 3] in the probability distribution for the 
mass ratio q = m2/mi [Eq. The remaining parame- 
ters are fixed and as described for our reference model in 
BKB. Most notably, we use a binary fraction /& = 1 (i.e. 
all stars are binaries) and solar metallicity {Z = 0.02). 

Synthesis Runs to Chosen Relative Accuracy: Since 
we explore an exceptionally large database of models to 
derive a fit for the inspiral rate, we need to minimize 
the considerable computational cost associated with each 
model. We do so by fixing the relative statistical ac- 
curacy of the inspiral event rate calculation from each 
model. 

Since the inspiral rates are linearly proportional to the 
number of relevant events n occurring in a given run, the 
relative accuracy of the inferred rates is proportional to 
^/n. For example, to obtain rates accurate within ^ 30%, 

^ To improve our computational efficiency, we restrict the Monte- 
Carlo generations to mi > 4 and miq > 4, since we are interested 
in just NS and BH in our present study. 



we require ~ 10. Therefore, for every model, we evolve 
each binary in succession, and we stop generating new 
binaries when either (i) a fixed number n of events has 
occurred or (ii) the total number N of generated binaries 
grows above a chosen threshold A'max ■ Our choices for 
the two specific thresholds, n and iVmax 7 vary depend- 
ing on the goal and nature of our calculation and are 
described in what follows. 

2.1. Example: Sampling the NS-NS Rate 

As an example of the above process, we choose random 
combinations of our 7 population synthesis model param- 
eters; for each combination, we select binaries according 
to the chosen initial distributions; we evolve each binary 
in succession and record the result at the end of the evo- 
lutionary sequence; we stop when we reach either n = 10 
NS-NS binaries which merge within the simulation time 
T — 10^" yrs^ or when we have sampled N = 10^ bina- 
ries, whichever comes first. Over the course of a month, 
on ~ 10 — 15 dedicated (of currently top-level speed) 
CPUs, we were able to evaluate 488 different population 
synthesis models. 

Each of these 488 runs provides a single, relatively low- 
accuracy (30%) but reasonable estimate of the NS-NS 
merger rate appropriate to the assumptions used in that 
run. This collection provides an unbiased sample of the 
NS-NS merger rate as a function of population synthesis 
model parameters in a 7-dimensional space. In Sec. 01 
we use these runs (along with other data) to generate 
fits, with which we can crudely and quickly estimate the 
NS-NS merger rate as a function of model parameters. 

3. ACCELERATING POPULATION SYNTHESIS 
SIMULATIONS 

While we successfully estimate the NS-NS merger rate 
with the direct approach outlined above for a large 
enough number of models, the case of the BH-BH 
merger rate is significantly more challenging because 
they are typically smaller than the NS-NS merger rates 
[cf. Fig. Q]. Conversely this implies that with the 
progenitor-generation scheme in the synthesis simula- 
tions the initial binary population is greatly dominated 
by systems which do not evolve to BH-BH binaries (and 
especially merging BH-BH binaries). 

We propose here that it is actually possible to use the 
above characteristics to our benefit and derive appropri- 
ate partitions (i.e., constraints) that eliminate the ma- 
jority of irrelevant progenitors, but at the same time do 
not eliminate a significant fraction of relevant progeni- 
tor. Although we concentrate here on the case of merging 
BH-BH binaries, it is clear that the method we describe 
below can be applied to any type of compact object bina- 
ries with low formation rates (relative to other binaries) . 

3.1. Partitions on the Initial Binary Parameter Space 
for Merging BH-BH Formation 

To accelerate the synthesis simulations for BH-BH 
merger rates we look generally for partitions in the parent 
parameter space: surfaces in the parameter space of all 
possible progenitors characterized by P = (mi, q, A, e) 

^ For brevity, we denote systems which merge due to gravita- 
tional radiation before the end of the simulated interval by "merg- 
ing binaries" . 
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that separate interesting from uninteresting initial bina- 
ries. 

To derive a partition well-suited towards finding pro- 
genitors of merging BH-BH binaries, we search for some 
combination of parameters C = (ci, C2, C3, C4, C5) such 
that the function 

p{P,C)=Ci logio(TOi) +C2 logio(m2) 

+C3logio(^) +C4e-C5 (3) 

is both (i) positive {p > 0) for all BH-BH binaries and 
(ii) negative {p < 0) for as many other non-BH-BH bi- 
naries as possible. [The form we use in Eq. (PJ for p is 
an arbitrary choice, having no motivation besides conve- 
nience.] Assuming such a partition can be found, we can 
safely ignore any progenitor parameters for which p < 0; 
doing so speeds up each run by a factor equal to the ratio 
of "weighted volumes" between the regions with p > 
and with p < 0.^ Clearly this process is beneficial only 
if the speed-up factor defined here is significantly higher 
than unity. 

Conceptually, these partitions can be found by a 
straightforward calculation. For example, we can take 
the recorded data from our 488 simulations and split the 
records into three groups: A (the set of all progenitor bi- 
naries which end up as merging BH-BH binaries), B (the 
set of all binaries which end up as merging NS-NS bina- 
ries), and E (everything else). Then we use some search 
algorithm to find those parameters C which have p > 
for as many members of A as possible and, furthermore, 
keep p < on a significant fraction of all members of B. 

To be quantitative, we use a robust, genetic-algorithm 
based search to find combinations C which maximize the 
following quantity: 

SiC) ^ ^ E C)) + ]^ E (^(-P(P' C-)) (4) 

^ ceA ^ ceB 

where 9{x) is a step function: 9{x) = 1 if x > and 
otherwise; and Na and Nb are the numbers of bina- 
ries in A and B, respectively. For example, the second 
term J2ceA^i~Pi^T^))/-^^ represents the fraction of 
elements of A which are misclassified by p. The ratio of 
the two prefactors, 100:1, is chosen so that a 1% error in 
classifying binaries of type A would count as significant 
as a 100% error in classifying binaries of type B. 

3.2. Requirements and Consistency tests 

The surprising fact is not that a partition like the one 
described above exists (i.e., not that a maximum of Eq. 
(@J can be found), but that it turns out to be both (i) 
useful: the ratio of the weighted volumes with p < 
and p > is surprisingly large, about 9:1, and so is the 
speed-up factor; and (ii) accurate: the fraction of binaries 
with P G A that are misclassified and their evolution is 
mistakenly ignored (i.e., p < 0) is very small (for our 
case, it is only 0.14%!). 

Of these two factors, the high accuracy is the most 
surprising. The partition we derive remains accurate 

^ Here the volume is weighted by the assumed probabihty distri- 
butions of the initial parameters, rather than geometric coordinate 
volume. 



TABLE 1 
Population Synthesis Runs 



Target Event 


Attempted 


Complete 


n 


-^max 


NS-NS(a) 


488 


418 


10 


10^ 


BH-BH 


312 


273 


10 


5 X 10"^ 


NS-NS{b) 


151 


123 


100 


5 X 10** 



for each model examined: we misclassify and ignore (i.e. 
p < 0) no more than a small fraction of the merging 
BH-BH binaries in any given run. Quantitatively, if we 
consider any simulation k in our set of 488 runs and de- 
note by Ak the set of all BH-BH binaries in that run, 
then J2ceAk ^{~P) < 0A2NAk', further, only 5% of runs 
have errors > 5%. 

3.3. Application: BH-BH runs via partitions 

As described above, we found a surprisingly robust par- 
tition which efficiently and accurately rejects the progen- 
itors of binaries which do not evolve into merging BH- 
BH binaries [cf. Eq. (j^J, using the coefficients given in 
the first line of Table [2]. Thus, to reduce the compu- 
tational burden needed to perform population synthesis 
runs geared towards the formation merging BH-BH bina- 
ries, we use this partition to augment the general proce- 
dure outlined in Sec. ??: we select progenitor parameters 
according to the BKB distributions; we reject some pro- 
genitors on the basis of this partition; we evolve each 
progenitor parameter combination in succession, stop- 
ping when we reach either n = 10 merging BH-BH bi- 
naries or = 5 X 10^ total binaries sampled, whichever 
came first. The second line in Table summarizes these 
choices. 

As in the NS-NS case, we chose random combinations 
of our 7 population synthesis model parameters. We 
rather quickly were able to evaluate 312 such combina- 
tions. 

3.4. Application: NS-NS runs via partitions 

As with BH-BH mergers, we can similarly attempt to 
employ partitions to reject systems which cannot possi- 
bly evolve into merging NS-NS binaries. Unfortunately, 
in sharp contrast to the BH-BH case, for NS-NS merg- 
ers we have not been able to dramatically accelerate 
each run - even with the use of several partitions si- 
multaneously (i.e., we reject a system if any partition is 
negative), each designed to filter out a specific type of 
contaminant (e.g., WD- WD binaries, disrupted binaries, 
WD-NS binaries, BH-BH binaries, etc; see Table |21 for 
the specific partitions used). Nonetheless, the identified 
partitions remained quite accurate (average error prob- 
ability 1.3%) and offered a non- negligible improvement 
in computation speed: a speed factor of 2.5, based on 
the ratio of "weighted" volumes] , we performed an addi- 
tional set of 151 runs, once again randomly distributed 
over model parameter space, designed to provide higher- 
accuracy (n = 100, leading to a statistical accuracy of 
10%) estimates of the NS merger rate. The third line in 
Table n] summarizes these choices. 

4. RESULTS: MERGER RATES AND MERGER RATE 
DISTRIBUTIONS 
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TABLE 2 
Partition table 



Target Type 


cl 


c2 


c3 


c4 


c5 


BH-BH 


0.971 


0.246 


-0.0167 


0.1567 


1.567 


NS-NS(b) 


+0.915 


+0.938 


-0.0915 


+0.1705 


+1.436 




+0.906 


+0.689 


-0.0415 


+0.1323 


+1.323 




-0.667 


+0.974 


-0.0293 


-0.0078 


-0.078 




-0.659 


+0.018 


-0.0293 


-0.1101 


-1.101 
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Fig. 1. — A smoothed, normalized histogram of log]^Q(7^), for TZ 
the merger rate for BH-BH (solid) and NS-NS (dashed) binaries, 
per year per galaxy similar to the Milky Way. The histogram has 
been smoothed over a length of 0.1 in the log, and normalized 
so J p(7?.)d logj^Q 7?. = 1. Our simulations are adequate to resolve 
the tail of each of these distributions. For example, for the BH-BH 
simulations (cf. Table 1), our simulations have adequate resolution 
to discover models with TZ > 10~^ yr~^. 



Tablensummarizes the relevant information for the ac- 
celerated runs we performed: the target event type (NS- 
NS or BH-BH), the number of runs (i.e., the number of 
parameter combinations used), the number of runs that 
were "completed" (i.e. which found the desired number 
of events of the target type),"* the number of target events 
n needed for a run to be considered "complete" , and the 
maximum number of primordial binaries allowed in these 
runs. Our sample consists of an order of magnitude more 
points than have been sampled before, distributed ran- 
domly throughout many plausible combinations of the 
seven model parameters varied. Each model is appro- 
priately consistent with observations of the Milky Way 
disk. 

Merger rate distributions: Figure ^ is a histogram of 
our results for the BH-BH and NS-NS merger rates from 
the complete sets of runs. In practice the two histograms 
represent a pair of probability distributions of BH-BH 
and NS-NS inspiral rates, assuming flat prior probability 
distributions for the population synthesis model param- 
eters. 

We can trust that these histograms accurately repre- 

* Almost all "incomplete" runs failed to find the maximum num- 
ber of merging binaries not because the model parameters would 
not produce them, but rather because practical matters (e.g. com- 
puter problems) prevented the run from finishing according to our 
desired conditions. 



sent probability distributions because we can accurately 
fit both rate functions over the 7-dimensional space of 
model parameters using the same data set. 

Seven-dimensional Fits for Merger Rates: Given the 
results on the NS-NS and BH-BH merger rates, sampled 
randomly through the seven-dimensional space of popu- 
lation synthesis model parameters, we attempt to obtain 
multi-dimensional fits. Such fits for the first time allow 
us to rapidly estimate these two rates without the need to 
resort to full population synthesis calculations that are 
typically highly demanding in computing power. The 
computational cost is dependent on the sophistication 
level of a given code. For StarTrack, a simulation for 10^ 
initial binaries with primary masses in excess of 4 M© re- 
quires about 80 hours on a single, AMD Athlon processor 
of top current speed. 

We have used several different fitting techniques (e.g., 
global polynomial fits at second and third order; non- 
parameteric fits using local quadratic approximations 
over the nearest 60 points; . . . ). All give comparable re- 
sults: based on their residuals, the derived fits turn out 
to be accurate to 50% (for NS-NS binaries) and 40% (for 
BH-BH binaries) [cf. Figure E]- In both cases, we find 
fits that have errors very close to the underlying errors 
in the data (due to statistical fluctuations in the small 
number n = 10 of merger events we observe),^ namely 
w w 30%. 

5. DISCUSSION 

In this paper we have systematically explored the de- 
pendence of predicted BH-BH and NS-NS merger rates 
on population synthesis model parameters by (1) esti- 
mating the rate for a number of different combinations of 
model parameters that exceeds earlier parameter studies 
by more than an order of magnitude; and then (2) fitting 
to the resulting data set, which provides us with a fast 
and moderately accurate (~ 0(40 — 50%)) estimate for 
these merger rates for general parameter combinations; 
the achieved accuracy is comparable to the underlying 
statistical accuracy of the runs used for the fit (30%). 

Separate from the focus of the present study, it is in- 
teresting to point out that our extensive database of 
synthesis models allows us to derive a prior probabil- 
ity distribution for the BH-BH merger rate, which has 
not been constrained empirically. In particular, this dis- 
tribution clearly implies that the BH-BH merger rate is 
higher TZ > 10~^yr~^ per galaxy similar to the Milky 
Way. Such a lower limit to the rate would lead to event 
rates for advanced LIGO sensitivity of at least tens of 
events per year. 

Furthermore, the derived fits lead to new inspiral-rate 
map that can immensely extend the exploration of the 
relevant parameter space. We can use them to address 
new questions that require a thorough understanding of 
how the inspiral rates vary with parameters, such as the 
following: 

1. Constraints on inputs (parameters): We have used 

^ The second NS-NS sample, for which n = 100, does not have 
enough points by itself to independently produce an accurate rate 
estimate. We used these points to test the fit obtained using the 
low-quality results, and found the errorrs were small. We also 
added these higher-accuracy points to our lower-accuracy data to 
marginally improve our overall fit. We found no significant change 
when we added these points. 
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Fig. 2.— A plot of our fits to the BH-BH (top) and NS-NS 
(bottom) merger rates versus our original estimate of that rate. 
Specifically, the plots show logj^Q 7?-j^(a;fc) versus logj^gT^j;, where 
7?.fc is the numerical value of the rate at and log7?.j is our fit to 
that rate data. In both cases, Tijix) was made from a quadratic fit 
to the nearest 60 points. The shaded region denotes our estimate 
for the relative uncertainty in each data point [i.e. the shaded 
region is the area between the curves 7?.(1 it l/vTO]. 



up to the development of the methods presented here. 
We expect each additional constraint will substantially 
reduce the set of population synthesis models consistent 
with observed astronomical populations. For example, 
as we will present in a subsequent paper, by merely com- 
paring a single prediction of population synthesis cal- 
culations (i.e., the NS-NS merger rate) with the known 
populations, we can exclude of order half our model 
space: half of our a priori model parameters are in- 
consistent with the observed distribution of binary pul- 
sars. Since population synthesis calculations implicitly 
produce a vast number of predictions which can be com- 
pared against observed populations (e.g., the supernovae 
rate; the statistics of X-ray binaries; etc) we expect that 
a more systematic study of experimental constraints will 
produce very stringent constraints on physically reason- 
able population synthesis parameters. 

Last we note that our code can and should be im- 
proved. For example, we are considering moving away 
from pure monte carlo, instead using nonrandom events 
(e.g., progenitor choices and supernovae kicks) and 
weighted systems to better sample the relevant parame- 
ter spaces and improve convergence. We are also explor- 
ing more sophisticated methods to seperate the progen- 
itors of "irrelevant" events from the progenitors of each 
target species. 

This work is partially supported by a NSF Gravita- 
tional Physics grant PHYS-0121416, a David and Lucile 
Packard Foundation Fellowship in Science and Engineer- 
ing, and a Cottrell Scholar Award from the Research 
Corporation to VK. 



only weak constraints on model parameters: we 
have assumed all values in some interval to be 
equally likely. Experimental results - most notably 
for the supernovae kick distribution - can provide 
probability distributions characterizing the likeli- 
hood of each parameter. Combining these distribu- 
tions with the inspiral rate, we can better estimate 
the NS-NS and BH-BH rate probability distribu- 
tions. 

2. Constraints on outputs (rates): We can use the 
empirically derived probability distribution for the 
NS-NS inspiral rate and map it hack to the 
population synthesis results. Such mapping will al- 
low us to construct probability distributions on the 
space of model parameters consistent with this em- 
pirical NS-NS distribution. We can further extend 
this concept and apply many additional observa- 
tional constraints simultaneously (e.g., on certain 
types of supernova rates, or on the formation rate 
of other known binary types, the rates of which can 
be empirically constrained; cf. Kim et al. 2004) 
to better constrain population synthesis model pa- 
rameters and results. 



We intend to undertake the above studies as a follow- 
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APPENDIX 

A. ABSOLUTE NORMALIZATION OF SYNTHESIS RUNS 

The scale factor s is a ratio s ~ Ng/N^ff between the number of steUar systems in the Milky Way and the number 
of stellar systems we have effectively sampled to select our n merging compact binaries. 

Effective sample size: The effective sample size Nf,ff is the number of stellar systems needed, on average, to produce 
N stellar systems with mi > 4, if all systems are drawn from an IMF which extends from the hydrogen burning limit 
m = O.OSMq to to = 1500: 

/.150 

Neff ^ N/ J dm (/)(m) . (Al) 

We use a Kroupa IMF: (/)(to) oc m~i-^ if to e [0.08, 0.SJM©, oc to"2.2 jf ^ g [o.5, 1]Mq, and oc to"^ '^ if to > IMq. 

Estimating the number of stellar systems in the galaxy: We choose Ng so that Ng times the average mass (according 
to our IMFs for toi and q and the binary fraction fb) of each stellar system (mtot) is equal to the total mass which 
should be formed in stars over the T = lOGyr lifetime of the Milky Way, MT: 



MT _ MT 
{mtot) {mi) (1 + fb {q}) 



^9 = 7Z-^ = 7Z7T7^^r7:JT- (A2) 



For the rate M at which mass is born in stars, we use the empirical estimate M « 3.5MQ/yr fRana 1 991tlBlitz 1997t 
ITacev 1985^) . The average values of toi and q are found from the Kroupa IMF and Eq. , respectively. 
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